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Abstract. - We present an analytical approach to the out of equilibrium dynamics of a 
class of kinetic lattice gases under gravity. The location of the jamming transition, the critical 
exponents, and the scaling functions characterizing the relaxation processes are determined. In 
particular, we find that logarithmic compaction and simple aging are intimately related to the 
Vogel-Fulcher law, while power-law compaction and super-aging behavior occur in presence of 
a power-law diffusion. 



Granular materials set in rapid motion by vibration Q] exhibit features such as universal 
velocity distribution and "de Gennes narrowing" — a physical signature of the cage 
effect H — which in spite of the non-equilibrium nature of the stationary state, closely 
resemble those observed in simple liquids at thermal equilibrium Q. In the opposite, quasi- 
static flow limit, slow compaction phenomena appear [n. During compaction, the free volume 
available to grains decreases, and the mobility steeply falls to zero, hence aging phenomena 
are expected to occur as is confirmed in several numerical simulations It has 

been suggested that in this regime a granular material should resemble a highly viscous liquid 



or a glass 1 13 , and several approaches have been proposed to describe different aspects of the 
granular compaction dynamics. These are mainly based on Langevin |p^| and Fokker-Planck 
equation [l5| , fluctuating nonlinear hydrodynamics |l6| , and mode-coupling theory [l7| . 

Two effects are responsible for the unusual behavior of a compacting granular material. 
First, collisions between the particles are inelastic, and energy has to be constantly pumped 
into the system. Second, at high packing density, stcric hindrance, and the associated cage 
effect, play a crucial role very similar to the one observed in amorphous systems. In this letter 
we shall be concerned precisely with this second effect, which in the case of gentle shaking is 
the dominant one. A number of lattice-gas models have been introduced to study numerically 
the slow dynamics induced by this effect |l^-pO|. Our approach allows to characterize in a 
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precise way the density relaxation and the aging dynamics in terms of the particle mobility. We 
show the existence of two relaxation regimes (exponential and aging) separated by a jamming 
transition. Using scaling arguments we investigate the asymptotic long-time behavior of the 
density and of the mean-square displacement, and compare the results with the numerical 
solution of the diffusion equation. We also check our analytic results with the Monte Carlo 
simulation of a gravity-driven version of the Kob- Andersen model (KA) Jl8], |lj . In this 
model, non-interacting particles are allowed to move only if a local kinetic constraint on the 
occupation of the nearest neighboring sites is satisfied. 

The model. - Consider a box of height H in contact with a particle reservoir of density 
p R located at z — H . All lengths are measured in units of the particle diameter. A constant 
gravitational field of strength g acts in the —z direction, forcing particles from the reservoir 
to enter into the box. The dynamical evolution of the local particle density p(z, t) is governed 
by the continuity equation, dtp{z,i) + d z J(z,t) = 0, with current given by the Fick's law, 
J(z,t) = —T(p)d z fj,(z,t). T(p) is the Onsager mobility and (J.(z,t) = ^ is the local chemical 
potential. For simplicity we shall suppose that the only interaction between the grains is due 
to the hard core repulsion. The exact Helmholtz free energy functional |23| for the lattice gas 
version of this model is 

F[p(z,t)} = k B T [ dz[{\- p)\n{\- p)+ p\np + 1 zp\ , (1) 
Jo 

where 7 = mg/k B T is the inverse gravitational length. In the case of granular materials the 
thermal energy of the grains is negligible so that T is a function of the externally imposed 
vibration intensity J27j] . It is neither the physical temperature nor the "granular temperature" 
usually associated with the average kinetic energy Q. For the KA model, the mobility T(p) 
vanishes as T(p) = Tq p (1 — p/p c )^, with the critical threshold p c = 0.88 for the simple cubic 
lattice jTsfl , and p c = 0.81 for the BCC lattice. In both cases = 3.1. A similar expression 
is found for the mobility of highly packed hard-sphere systems [^J, ^5|. One should note, 
however, that the free energy functional given by Eq. (|l|) is exact only for the lattice gas. 
In the case of hard spheres it can be replaced by an approximate expression derived from 
the Carnahan-Starling equation of state [p3|,p6[, without any qualitative modifications to the 
theory. Substituting T(p) into the continuity equation we are lead to, 



dp(z,t) d 
dt dz 




1 dp 
-, jr +7P 

1 — p oz 



(2) 



where the time is now measured in units of l/T^k^T. We should stress that the local density 
approximation in which the functional form of mobility is assumed to remain valid in an 
inhomogeneous system is highly non-trivial. One of the goals of this paper is to assess to 
what extent this approximation works in the case of kinetic lattice-gas models. 

The boundary conditions require vanishing of the current at z = 0, J(0,t) = 0; and 
p{H, t) = p R . The stationary state is obtained when dp/dt = 0, which implies that J(z, 00) = 
for all z. The jamming transition corresponds to the locus in the parameter space (7, p R ), at 
which a layer of critical density first appears, p(0, 00) = p c . This happens when 



7c (p R ) = -^hi 



Pc(l - Pk) 



Pb (1 - Pc ) 



(3) 



Depending on the value of 7 two very different types of stationary profiles are found. In Fig. [I] 
we compare these with the stationary profiles found from the Monte Carlo simulation of the 
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gravity-driven KA model on the BCC lattice |2l| . A perfect agreement is obtained with no 
adjustable parameters. It is important to note that for 7 > j c (Pr) the stationary profiles 
are no longer equilibrium distributions, i.e. they do not minimize the Helmholtz free energy 
functional Eq. |]. 




Fig. 1 - Examples of the stationary profiles above and below the critical temperature (solid lines). 
For 7 < 7c (pr) (high temperature) the equilibrium density profile is Poo(z) = l/[exp(7z + ?7) + 1], with 
r\ = In [(1 — /5r)/pr] — 7-ff. At low temperatures, 7 > 7 c (pr), the stationary profile is poo(z) = p c for 
z < zo(pn) = H [1 — 7c(Pr)/t]; and p x {z) = l/[exp(7z + r/) + 1] for z > zo(Pr). The squares are the 
asymptotic densities obtained through simulation of the gravity-driven KA model on the BCC lattice 
(see pi| for details on simulation). 



High-temperature phase, 7 < 7 c (/0r). - At high temperature we expect an exponential 
relaxation towards the equilibrium distribution, p(z,t) x poo(z) + g(z)e~ t ^ T . We have solved 
Eq. (U) numerically and checked that the approach to equilibrium is indeed exponentially fast 
for any 7 < 7 c (pr). To see how the characteristic time r depends on the various parameters 
of the system we perform a scaling analysis of Eq. (||) . We first note that associated with the 
two right-hand terms of Eq. (^) there are two temporal scales. These can be identified as the 
diffusion time t± oc H 2 , and the drift time T2 oc H/7. The characteristic time must, therefore, 
satisfy an exact equation t _1 = jjjjJ-^H; p K ), where J-~(x; p R ) is a scaling function. The 
reason for the prefactor 7r 2 /4 will become clear from the forthcoming analysis. To further 
explore the properties of Eq. (Q) we now study its linear version, 

d P (z,t) _ d 2 P dp 

dt dz 2 1 dz' v ' 

For Eq. (^) , the characteristic time of approach to equilibrium can be related to the singular- 
ities of the temporal Laplace transform of the density profile. For r > 4/7 2 the characteristic 
time is found to satisfy 



exp W7 2 -4t-i} = 7 + v/ /' £=j , (5) 

while for r < 4/7 2 it satisfies 

cos [h^/At- 1 - 7 2 } = r 7 2 /2 - 1 . (6) 
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We note that for 7 = 0, r 1 = 7r 2 /4_ff 2 . For finite 7 the inverse relaxation time for linear 
Eq. can be written in a scaling form, t _1 = -^^f(jH). Unlike the scaling function 
T(x\ p R ) for the non-linear Eq. (||), f(x) does not depend on the reservoir density. For x = 0, 
/(0) = 1; for large x, f(x) x (4/7r 2 )a; 2 exp(— 2). The graphs of /(x) and of its asymptotic 
form are presented in Fig. |2|. 




Fig. 2 - Scaling function f(x) (solid curve) and its asymptotic form (dashed line) as discussed in text. 

Fig. 3 - Inverse relaxation time, r , for 7 = 0.1 and H = 10 with p c = 0.88 and <j> = 3.1, 
corresponding to the Kob-Andersen model on a simple cubic lattice. The points are the result of 
numerical integration of Eq. (Q). Inset shows the same data on the log scale. The characteristic time 
diverges with exponent <j> — 2 as the jamming transition is approached. 

With the insight gained from the study of Eq. (Q) we can now explore the full non-linear 
Eq. (|^). We first consider the special case of zero gravity pS[ ]. For 7 = particles diffuse freely 
from the reservoir until a uniform density profile Pac(z) = p R is established. At a low reservoir 
density p R <C p c , we expect that the non-linearities of Eq. (0) should be irrelevant and the 
relaxation time should reduce to t' 1 rs jjj?- Linearizing Eq.(g) around the equilibrium state 
we find that the relaxation time for 7 = is 

or equivalently ^"(0;p R ) = (1 — p R /p c )^/(l — p R ). Eq. (0) is in perfect agreement with the 
numerical integration of Eq. Q). As expected, the relaxation time diverges as p R — » p c . The 
exponent characterizing this divergence is <j). 

Finally in the presence of a gravitational field we find, by numerical integration of Eq. @ , 
that as 7 — > 7c(/0r), the density of the first layer approaches p c , poo(Q) Pc and T ~ 
(7c — l)^~ 2 ■ The relaxation time diverges with exponent (/> — 2, see Figj3L implying that the 
dynamics is faster than in the zero gravity case. Comparing with eq. (|7j), we see that the 
jamming transitions in the homogeneous and inhomogeneous systems, therefore, belong to 
distinct universality classes. 

Low-temperature dynamics, 7 > "f c (p R ). - At low temperatures, the density of the lower 
layers closely approaches the critical threshold, p c . At such large packing fractions the move- 
ment of particles is highly restricted, and even the slightest increase of density requires a 
rearrangement of a huge number of grains. In this regime we expect slow relaxation and aging 
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phenomena to appear. To check this supposition one should solve Eq. (Q). Unfortunately, due 
to its highly non-linear nature no analytic solution is possible. Nevertheless an asymptotic 
solution can be found. Let us focus our attention on the bottom layers, z < Zo(p R ), of the 
suspension, Fig. [l| At sufficiently long times, the density of these layers will be close to the 
critical one, p{z,t) ~ p c - To the lowest order in A(z,t) = 1 — p(z,t)/p c , Eq. (0) simplifies to, 



dA(z,t) 
dt 



dA* 



-7" 



dz 



To solve this non-linear equation we propose a scaling ansatz A(z, t) — A(zjt c 
into Eq. M) we see that this form is a solution if a = 1 and 



(8) 

Substituting 



A(z,t) 



</>7 1 



(9) 



To check the asymptotic solution, Eq. (g), we have numerically integrated Eq. (g) for 7 > 
Jc(Pr)- For large times a perfect agreement between the numerical solution and the Eq. (^J) 
is found, as can be seen from Fig. [|. In the absence of gravity, 7 = 0, relaxation is slower 
and characterized by a different dynamic exponent, A(z,t) ~ t^ 1 ^ |p8|| . The same behavior 
is observed in the KA model (see ref. |2j| and Fig. ||). 

Up to now our discussion has been motivated by the dynamics of the KA model which is 
characterized by a power law mobility coefficient. However, there are various systems for which 
the mobility vanishes according to the Vogel-Fulcher law, = Tq p exp [a p c /(p — p c )]. In 
this case an analysis similar to the one presented above leads to an asymptotic scaling solution, 
which for very large times reduces to 



A(z,t) = 



\n(t/z) 



(10) 



Interestingly, a similar logarithmic law was discovered in the granular compaction experi- 
ments p9| , and is supported by several analytical approaches p(i|-]34|] . 
We now turn to the discussion of the aging phenomena [MIS 



36 . Consider the two-times 



mean square displacement of grains, B(t,t w ). If t is sufficiently larger than i w , this can be 
written as 

B(t,t w ) = dsT[p(z,s)] . (11) 
For power-law diffusion and z < zo(p R ) we find, to leading order in t and t w , 

s(t,tw)~*ir M -t 1-M , (12) 

with the exponent p = <j)/{4> — 1). Usually <f> > 1, so that p > 1, which corresponds to a 
super-aging regime [Q. This means that the effective structural relaxation time grows as t£, 
what is faster than the age of the system, t w Qj. This scaling behavior has been observed 
in the simulation of the gravity-driven KA model |39[ . The comparison with the zero-gravity 
case, for which a simple aging (t/t w ) was found J28[ , once again shows the important role 
played by gravity. 

For the Vogel-Fulcher law we find, to the leading order in t and t w (and z < Zo(p R )), a 
simple aging scenario, B(t,t w ) ~ log (i/i w ). It is interesting to observe the different behavior 
of this function and JT^ ) at finite waiting times i w . In the former case, lim^oo B(t, t w ) = 00, 
i.e., a weak ergodicity scenario |37j]; while in the latter, a finite limit is obtained (which, 
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Fig. 4 - Power- law density relaxation for 7 > j c {pr)- The points are result of numerical integration 
of the Eq. (^), while the solid line is the asymptotic solution given by Eq. (^|). 

Fig. 5 - Density relaxation found in the simulation of the gravity driven KA model on the BCC 
lattice. Averages are over 100 samples and H = 4L — 80, pn — 0.1, and 7 = 0.1. The line is the 
theoretical prediction t -1 /W -1 ). 

however, vanishes as t w — > 00). Nevertheless, the manner in which time-translation invari- 
ance is violated is, in a sense, similar. Indeed, if we consider the "triangle relation" f38|| , 
B(ti,t^) — f \B(t\,t2), B{t2,H)], where the times tx, t2, and £3 are in increasing order, it is 
straightforward to check that in both Vogel-Fulcher and power law cases f(x, y) = x+y, imply- 
ing that displacements over non-overlapping time intervals are statistically independent. This 
feature does not hold in the presence of activated aging for which B(t,t w ) ~ log</logt w |8). 

Conclusions. - To summarize, we have presented an analytical study of the dynamics 
of a class of kinetic lattice models. Our approach based on the local density approximation 
for the Onsager mobility allows us to predict the density relaxation law and the nature of 
aging behavior. The results reproduce successfully the behavior observed in the numerical 
simulation of the KA and gravity-driven KA model, and they could be relevant to systems for 
which steric hindrance and cage effect are dominant, such as polydispersed colloidal suspension 
in gravitational field (2^,^||, and slowly compacting dense granular matter. For hard-sphere 
systems it is known that the diffusion coefficient vanishes as a power-law with </> between 2 
and 3, and the sedimentation profile is approached very slowly (one year, roughly) [Q; hence 
a super-aging behavior in the two-times mean-square height displacement with exponent /i 
between 3/2 and 2 should be observed. 

if if. if. 
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